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We use first principles calculations to study structural, vibrational and superconducting properties 
of H 2 S at pressures P > 200 GPa. The inclusion of zero point energy leads to two different possible 
dissociations of H 2 S, namely 3 H 2 S —* 2 H 3 S + S and 5 H 2 S —* 3 H 3 S + HS 2 , where both H 3 S and 
HS 2 are metallic. For H 3 S, we perform non-perturbative calculations of anharmonic effects within 
the self-consistent harmonic approximation and show that the harmonic approximation strongly 
overestimates the electron-phonon interaction (A « 2.64 at 200 GPa) and T c . Anharmonicity 
hardens H-S bond-stretching modes and softens H-S bond-bending modes. As a result, the electron- 
phonon coupling is suppressed by 30% (A « 1.84 at 200 GPa). Moreover, while at the harmonic level 
T c decreases with increasing pressure, the inclusion of anharmonicity leads to a T c that is almost 
independent of pressure. High pressure hydrogen sulfide is a strongly anharmonic superconductor. 


Cuprates [1] have for many years held the world 
record for the highest superconducting critical tempera¬ 
ture (T c = 133 K) [2]. However, despite almost 30 years 
of intensive research, the physical mechanism responsible 
for such a high T c is still elusive, although the general 
consensus is that it is highly non-conventional. The dis¬ 
covery by Drozdov et al. [3] of T c = 190 K in a diamond 
anvil cell loaded with hydrogen sulfide (H 2 S) and com¬ 
pressed to about 200 GPa breaks the cuprates record and 
overturns the conventional wisdom that such a high T c 
cannot be obtained via phonon-mediated pairing. 

The claim that hydrogen at high pressure could be 
superconducting is not new [4] and it was recently sup¬ 
ported by first principles calculations based on the har¬ 
monic approximation applied to dense hydrogen [5-8] 
and several hydrides [9-15]. More recently, two theo¬ 
retical papers predicted the occurrence of high T c su¬ 
perconductivity in high-pressure sulfur-hydrides [16, 17]. 
However, as shown in Refs. [18, 19], anharmonicity can 
be crucial in these systems. For example, in PdH, the 
electron-phonon coupling A parameter is found to be 1.55 
at the harmonic level, while a proper inclusion of anhar¬ 
monic effects leads to A = 0.40 [18], in better agreement 
with experiments. Thus, in hydrogen-based compounds, 
the phonon spectra are strongly affected by anharmonic 
effects. 

Several first principles calculations [16, 17, 20, 26] sug¬ 
gested that decomposition of the H 2 S sample occurs 


within the diamond-anvil cell at high pressures. The 
high-T c superconducting material is therefore very un¬ 
likely to be H 2 S, while H 3 S is the obvious candidate for 
the H-rich decomposition product. 

Here we study the structural, vibrational and super¬ 
conducting properties of H 2 S above 200 GPa, where the 
highest T c occurs. We show that the inclusion of zero 
point motion in the convex hull at 200 and 250 GPa sta¬ 
bilizes two metallic structures, H 3 S and HS 2 . Finally, 
we show that, contrary to suggestions in previous work 
[16, 20], the harmonic approximation does not explain 
the measured T c in H 3 S, and the inclusion of anharmonic 
effects is crucial. 

As decomposition of H 2 S has been demonstrated in the 
experiments of Ref. [3], it is crucial to develop an under¬ 
standing of the different H/S compounds that might be 
stable in the pressure range of interest. We therefore per¬ 
form a search over 43 H/S stoichiometries, determining 
the stoichiometries at which stable structures exist, and 
the associated crystal structures. These searches were 
performed using the ab initio random structure search¬ 
ing (AIRSS) method [21, 22] and the CALYPSO particle 
swarm optimization method of Ref. [23]. More informa¬ 
tion about the searches is provided in the Supplemental 
Material [24]. 

The results of the structure searching are shown in Fig. 
1. At 200 GPa, without zero-point energy (ZPE), the 
only energetically allowed decomposition is 3H 2 S —► 2H 3 S 



2 



Composition (S fraction) 

FIG. 1. Results of structure searching at 200 and 250 GPa. 
Convex hulls are shown as continuous lines, with and without 
the inclusion of zero point energy (ZPE). 


+ S, in agreement with previous calculations [16, 20, 25, 
26]. H 3 S crystallizes in the space group Im3m, as shown 
in [16]. When ZPE is included, a second decomposition 
becomes possible at 200 GPa, namely 5H 2 S —> 3H 3 S + 
HS 2 , where HS 2 crystallizes in a structure of space group 
C2/c with 12 atoms/cell. At 250 GPa and above, the 
latter decomposition is allowed even without ZPE, and 
the C2/c HS 2 structure undergoes a phase transition to 
a more stable C2/m structure with 6 atoms/cell. Each 
of the HS 2 structures is metallic. Finally, at 300 GPa, 
an HS phase becomes stable [24]. Detailed information 
on the crystal structures is provided in the Supplemental 
Material [24]. 

Having determined the most stable crystal structures 
at high pressure, we turn to the study of vibrational prop¬ 
erties [27, 28]. We consider the Im3m H 3 S structure at 
200 GPa [33]. In this structure each H atom is twofold 
coordinated and has 6 neighbors, 2 of which are S atoms 
while the other 4 are H atoms. H vibrations can then 
be decomposed into HS bond-stretching modes (Hy), in 
which an H atom moves towards one of the two S atoms, 
and bond-bending modes (HjJ, in which one H atom 
moves in the direction perpendicular to the H-S bond 
(see Fig. 1 in [24]). The harmonic phonon spectrum of 
H 3 S is shown in Fig. 2 and overall shows a clear sepa¬ 
ration into H modes at high energy and S modes below 
75 meV. To gain more insight we use Wannier interpola¬ 
tion [31, 32] of the electron-phonon matrix elements and 
evaluate the electron-phonon contribution to the phonon 



FIG. 2. Phonon dispersion, phonon density of states projected 
onto selected atoms and directions, and the Eliashberg func¬ 
tion of H 3 S in the harmonic approximation (top) and with the 
inclusion of anharmonic effects (bottom) for H 3 S at 200 GPa. 
H_l and Hy label displacements of an H atom in the directions 
perpendicular or parallel to a H-S bond. The magnitude of 
the phonon linewidth is indicated by the size of the red error 
bars. 

linewidth, as [34]: 

lffnm( k ^ k + q)| 2 ^(ekn)^(£k+qm) • (1) 

k k,n,m 

Here w q „ are the phonon frequencies, N k the number of 
electron-momentum points in the grid, g£ m (k, k + q) = 
(kn|<5Vftrs/<fot ql ,|k + qro) is the electron-phonon matrix 
element, Vks is the Kohn-Sham potential, and u q!y is a 
phonon displacement. The Kohn-Sham energy and eigen¬ 
functions are labeled £kn and |kn). The electron-phonon 
coupling at a given phonon-momentum q for a phonon 
mode v can be obtained [34] from the phonon linewidth 
as A q „ = 2 nJ?N(o) • 
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As shown in Fig. 2, at the harmonic level, the phonon 
linewidths of the H vibrations is fairly uniform through¬ 
out the spectrum. The contribution of each mode to the 
average electron-phonon interaction, A = A q v/N q , 

can be obtained from the isotropic Eliashberg function 


a 2 F{u>) 


1 

2 W q 


A q ^Wq^l5(w 
q v 


Uqy) 


( 2 ) 


where N q is the number of phonon-momentum points in 

the grid. A(w) = 2 a ^ dto' and then A = A(oo). 
We find A = 2.64 (see Table IV), which is larger than 
that obtained in Refs. [16, 20, 26] with a much coarser 
sampling of the BZ. This huge value of A comprises sub¬ 
stantial contributions from many H vibrational modes. 
The situation is therefore very different from MgB 2 in 
which a single mode dominates A. 

Given the low mass of H and the consequent large 
phonon displacements, we investigate the occurrence of 
anharmonic effects using the stochastic self-consistent 
harmonic approximation (SSCHA) developed by some of 
us [18, 19, 35]. As shown in Fig. 2 (bottom), the an¬ 
harmonic correction leads to non-trivial changes in the 
harmonic spectrum. While it is very clear that all H 
bond-stretching modes are hardened, the effect on H 
bond-bending modes is less straightforward. By com¬ 
puting the average phonon frequency as of Hy and of Hj_ 
modes we find that u7y ar « 158.1 meV and Wy nh « 203.3 
meV, while for bond-bending modes w l [ ar « 157.0 meV 
and w anh « 147.9 meV. Thus bond-stretching modes are 
hardened, while bond-bending modes are softened. 

It is important to remark that the large and most 
dispersive mode along PT is strongly hardened at the 
anharmonic level and undergoes a non-trivial change in 
polarization, as can be seen from the large effect of an- 
harmonicity on the phonon linewidth 7 q „ in Fig. 2. It 
is worhwhile to recall that the phonon-linewidth depends 
on the phonon eigenvector but not on the phonon energy. 
This effect demonstrates the need to calculate not only 
the phonon frequencies at the anharmonic level, but also 
the phonon polarizations. 

The anharmonic electron-phonon interaction is A = 
1.84, which is 30% smaller than the harmonic result. This 
reduction is mostly explained by the hardening of the Hy 
modes. In contrast to the harmonic case which shows 
uniform coupling over all modes, the anharmonic Eliash¬ 
berg function has two main peaks, a broad peak in the 
40-75 meV region, and a second one in the 175-200 meV 
region. Their contributions to A are 0.59 and 0.77, re¬ 
spectively, accounting for « 73% of the total A. We note, 
however, that the logarithmic average of the phonon fre¬ 
quencies, wiog, is only weakly enhanced by anharmonicity 
(see Table IV). 

The superconducting critical temperature can be ob¬ 
tained either from the McMillan equation or the isotropic 
Migdal-Eliasliberg approach. However, it is well known 


[36] that the use of the McMillan equation for such val¬ 
ues of A leads to a substantial underestimation of T c . We 
solved the isotropic Eliashberg equations [24] and found, 
contrary to claims in previous publications [16, 20 ], that 
calculations based on the harmonic phonon spectrum do 
not explain the measured T c as, even using large val¬ 
ues of n* [37, 38], T c is substantially overestimated (i.e., 
T c = 250 K for /r* = 0.16 [24]). When the anharmonic 
phonon spectrum and electron-phonon coupling are used, 
the Migdal-Eliashberg equations account for the experi¬ 
mental T c when the value /i* = 0.16 is used, as shown in 
Table IV. The superconducting gap at zero temperature 
is A « 36.5 meV. 

Interestingly, the large anharmonic effects lead to very 
different variation of T c with pressure. By repeating the 
calculation for the Im3m structure at 250 GPa, we found 
at the harmonic level and using the Migdal-Eliashberg 
equations with the same values of /j,* = 0.16, that T c = 
226 K, decreasing with increasing pressure. However, at 
the anharmonic level we find T c = 190 K, essentially 
independent of pressure in the region 200-250 GPa. 

Finally, we consider the extent to which the occur¬ 
rence of large anharmonic effects can explain the isotope 
shift in D 2 S. At 164 GPa, T C (D 2 S) = 90 K, leading to 
an isotope coefficient a ss 1.07, which is substantially 
enhanced from the canonical BCS value of a ss 0.5. As¬ 
suming a similar decomposition of D 2 S into D 3 S and S 
at high pressures, we calculate the anharmonic phonon 
spectrum (see [24]) and electron-phonon coupling in D 3 S 
at 200 GPa. We find at the anharmonic level that the 
electron-phonon coupling is essentially unaffected, while 
wi 0 g is softened from 92.9 meV to 73.3 meV, leading to an 
isotope coefficient of a = 0.35, which is strongly reduced 
from the BCS value but inconsistent with the value of 
a « 1.07 found in experiments. Thus, contrary of what 
claimed in Ref. 39, anharmonicity reduces a. 


TABLE I. Electron-phonon interaction and logarithmic av¬ 
erages of phonon frequencies, with and without anharmonic 
effects. The T c s are calculated using the isotropic Migdal- 
Eliashberg equations (T^ E ). A value of fi* = 0.16 is used. 
Data for T c calculated with the McMillan equation is pro¬ 
vided in the Supplemental Material [24]. Frequencies are in 
meV and T c s are in K. 


Compound 

A har. 

, ,har 
W lo E 

A a„h 

"Ef 

rjiMtl,har 

rpMtl,anli 

c 

T c (Exp) 

H 3 S (200 GPa) 

2.64 

90.4 

1.84 

92.86 

250 

194.0 

190 

H 3 S (250 GPa) 

1.96 

109.1 

1.71 

101.3 

226 

190 


D 3 S (200 GPa) 

2.64 

68.5 

1.87 

73.3 

183 

152.0 

90 


We have studied the structural, vibrational and super¬ 
conducting properties of high pressure H 3 S. We have in¬ 
cluded zero point motion when comparing the stabilities 
of different H/S phases, which has been neglected in other 
publications so far. This is important because zero point 
motion stabilises a new phase at P > 200 GPa. In par¬ 
ticular, we found that decomposition of HS 2 into metallic 
phases can occur following two main paths, namely 3H 2 S 
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—> 2 H 3 S + S and 5 H 2 S —> 3 H 3 S + HS 2 . We have per¬ 
formed a detailed study of the vibrational properties of 
high pressure H 3 S and D 3 S, finding that the phonon spec¬ 
tra are strongly affected by anharmonic effects. Anhar- 
monicity hardens H-S bond-stretching modes and soft¬ 
ens H-S bond-bending modes. Moreover, anharmonicity 
leads to a reduction in the electron-phonon coupling by 
ss 30% and to an approximately constant T c in the range 
200-250 GPa. Our work demonstrates that the super¬ 
conducting properties of high pressure H 3 S can only be 
properly described by including both nuclear quantum 
effects and anharmonicity. 
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Technical details for the structural searching 


We use the CALYPSO and AIRSS codes. Structure searching were performed at 200, 250, and 
300 GPa for several H-S compounds (H8S, H7S, H6S, H5S, H4S, H3S, H2S, HS, H3S2, H2S3, H5S2, 
H9S2, HS2, H4S3, HS3, H3S2, H7S2) with maximum eight formula unit in the models. Each 
generation contained 40 structures, and the first generation was produced randomly with symmetry 
constraint. All structures were locally optimized using density functional theory with the Perdew- 
Burke-Ernzerhof (PBE) [1], generalized gradient approximation implemented in the Vienna ab initio 
simulation package[2, 3]. An energy cutoff of 700 eV and a Monkhorst-Pack Brillouin zone sampling 
grid with a resolution of 0.5 A ” 1 were used in structure searches. The 60% lowest-enthalpy structures 
of each generation were used to produce the structures in the next generation by local PSO technique, 
and the remaining 40% structures were randomly generated within symmetry constraint to enhance 
the structural diversity. Typically, the structure searching simulation for each composition was 
stopped when 1000 successive structures were generated after a lowest energy structure was found. 
A number of distinct low-enthalpy structures found were then re-optimized with denser grids better 
than 0.2 A -1 and a higher energy cutoff of 1000 eV. The lowest-enthalpy structures were then chose 
to draw the convex-hull. 
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Convex hull at 300 GPa 



Composition (S fraction) 

FIG. 3. Results of structural searches at 300 GPa. The continuous line shows the convex hull. 


Effect of the zero point energy on the DtS convex hull at 200 GPa 



Composition (S fraction) 


FIG. 4. Results of structural searches at 200 and 250 GPa for D-S structures including zero point energy. The continuous lines 
show the convex hull 





Crystal structures for H and S 


The crystal structure for Hydrogen is from Ref. [5]. The crystal structure for S is from Ref. [6] 


Crystal structures of HS 2 at 200, 250 and 300 GPa 


TABLE II. Crystallographic data for HS2 at 200, 250 GPa as obtained from structural searches. At 250 and 300 GPa, HS2 
adopts the same C2/m structure as at 250 GPa. 


Structure Pressure 
(GPa) 

Structural parameters 
(A, deg.) 

Atomic Positions 


C2/c 

200 

a=6.7827, b=4.1876, c=7.5401 S 8f 

0.07702 

0.12708 0.42711 



a = 90, /? = 137.7464, 7 = 90 S 8f 

0.89224 

0.37399 0.79574 



H 8f 

0.28345 

0.86882 0.42463 

C2/m 

250 

a=7.2073, b=2.947, c=3.6324, S 4i 

0.33820 

0.5 0.51566 



a = 90, P = 60.2287,7 = 90 S 4i 

0.59068 

0.0 0.14405 



H 4i 

0.61860 

0.5 0.94147 


The HS 2 crystal structures found with the CALYPSO and AIRSS codes are shown in Table III. 
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Crystal structures of HS and at 200, 250 and 300 GPa 


TABLE III. Crystallographic data for HS at 200 and 300 GPa as obtained from structural searches. At 250 GPa and 300 GPa, 
HS 2 adopts the same C2/m structure. 


Structure Pressure Structural parameters 



Atomic Positions 



I4i /amd 

200 a=b=2.9399, c=9.0531 

S 

16h 

0.0 

0.0 

0.27894 



H 

I 6 I 1 

0.5 

0.0 

0.62786 

C2/m 

300 a=9.4579, b=2.7388, c=2.749, 

S 

4i 

0.41556 

0 

0.83415 


a = 90 , /3 = 73.1325, 7 = 90 

S 

4i 

0.33437 

0.5 

0.41585 



H 

4i 

0.00217 

0.0 

0.24756 



H 

4i 

0.25205 

0.5 

0.99709 


Crystal structure of H 3 S at 200 and 250 GPa 



FIG. 5. Crystal structure of H 3 S at 200 and 250 GPa. Hydrogen is depicted in cyan while sulfur is in yellow. The red (green) 
arrows label the H-S bond-bending (-stretching) modes, labeled Hx (H||). The volumes at 200 and 250 are 13.3334 A 3 , and 
12.4925 A 3 , respectively. 


Electronic structure of H 3 S at 200 GPa 


The electronic structure of H 3 S at 200 GPa is shown in Fig. 7 in the fat-bands representation 
including decomposition into H and S atomic states. In a 20 eV energy window around the Fermi 
level (e_f), the electronic structure can be fairly well interpreted in terms of free electrons on a bcc 
lattice. However, in the proximity of the Fermi level, there is a substantial hybridization between the 
H and S electronic states, leading to avoided crossings at special points N, T and along the H-N high 
symmetry direction. The hybridization results in a peak in the density of states at approximately 
0.17 eV below e/. At precisely e/, the density of states per spin is N( 0) = 0.33 states/eV/spin/H 3 S 
cell, which is essentially identical to the free-electron value. 










10 



FIG. 6 . Left: Free-electron bandstructure and density of states calculated for a bcc lattice with the same lattice parameter 
as H 3 S at 200 GPa. The electronic structure and density of states of H 3 S at 200 GPa are also shown for comparison. Right: 
Electronic structure of H 3 S. The thickness of the band is proportional to the projection of the electronic state over a chosen 
atomic orbital (fat bands representation). 



FIG. 7. Fermi surface of H 3 S at 200 GPa including 5 sheets (right). The hole pockets at T are visible in the center. The large 
Fermi surface contributing most of the density of states is shown in transparent colors. 


The Fermi surface is composed of 5 sheets. The avoided crossing at the zone center generates three 
hole pockets centered at T and a cubical electron Fermi surface centered at the N point. Finally, an 
additional large Fermi surface sheet arises from the electron-pocket at H. 


Vibrational and superconducting properties of D 3 S at 200 GPa 


The phonon dispersion of D 3 S calculated using harmonic linear response theory and the stochastic 
self-consistent harmonic approximation is shown in Fig. 8 . 
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FIG. 8 . Phonon spectrum, Eliashberg function and integrated electron-phonon coupling of D 3 S at the harmonic and anharmonic 
levels (SSCHA). 


Vibrational and superconducting properties of H 3 S at 250 GPa 


The phonon dispersion of D 3 S calculated using harmonic linear response theory and the stochastic 
self-consistent harmonic approximation are shown in Fig. 9. 
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FIG. 9. Phonon spectrum, phonon density of states projected onto selected vibrations and the Eliashberg function and 
integrated electron-phonon coupling of H 3 S at 250 GPa at the harmonic and anharmonic levels (SSCHA). 
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Migdal-Eliashberg 


We solve the Isotropic Migdal-Eliashberg (ME) equations using either the harmonic Eliashberg 
function or the one calculated within the SSCHA. The equations are solved in the Matsubara fre¬ 
quency space using 512 Matsubara frequencies. The superconducting gap is obtained from the lowest 
Matsubara gap A n=0 = A and is plotted in Fig. 10. 



FIG. 10. Superconducting gap (A) as a function of temperature from the solution of the isotropic Migdal-Eliashberg equations 
applied to H 3 S. The continuous lines refer to anharmonic phonon (SSCHA) while the harmonic phonon are denoted by dashed 
lines. 


Superconducting properties using different approximations 


TABLE IV. Electron-phonon interaction and a logarithmic average of phonon frequencies using the SSCHA, and without 
anharmonic effects. The T c s are calculated using the SSCHA phonon spectrum and with either the McMillan equation (T^ M ) 
or by solving the isotropic Migdal-Eliashberg equations (T^ E ). A value of /r* = 0.16 is used. 

Compound A har tgg (mcV) A anh (mcV) 7^ 1M ' hal ' 7) MM - anl1 T^ IE - har T^ E - anh T c (Exp) ~ 


H 3 S (200 GPa) 2.64 

90.4 

1.84 

92.86 

158.8 

124.9 

250 

194.0 

H 3 S (250 GPa) 1.96 

109.1 

1.71 

101.3 

155.25 

127.2 

226 

190 

D 3 S (200 GPa) 2.64 

68.5 

1.87 

73.3 

120.4 

100.3 

183 

152.0 


90 
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Effects of the vibrational energy on the pressure 


In Fig. 11 the contribution of the atomic vibrations to the total energy is shown in the quasi¬ 
harmonic approximation and SSCHA. The ground state total energy, without the zero point en¬ 
ergy (ZPE), is fitted to a second order polynomial as a function of the volume V, E 0 (V) = 
Aq-\-BqV + CqV 2 - The fitting parameters are reported in Table V. We add the vibrational contribu¬ 
tion to Eq(V) calculated in the quasiharmonic approximation and in the SSCHA. This contribution 
is calculated at two volumes, those shown in the right panel of Fig. 11. The vibrational energy E v 
is fitted linearly to E v (V) = A v + B v V. As shown in Fig. 12, the linear fit provides a very good 
approximation to E v (V). The linear form of E v (V) is obtained within the quasiharmonic approxi¬ 
mation for H 3 S and D 3 S, and in the SSCHA for H 3 S. The ZPE energy is obtained using a 6x6x6 
phonon mesh. 

The correction to the pressure from including the atomic vibrations is calculated from the E(V) = 
E 0 (V) + E v (V) curves and P(V) = —&E(V)/AV. These curves are shown in Fig. 11. Neglecting 
the ZPE underestimates the pressure. The correction to the pressure from the ZPE is smaller for 
D 3 S than for H 3 S due to the smaller ZPE of deuterium. The SSCHA gives a small correction to the 
pressure obtained within the quasiharmonic approximation. The pressure corrections to the volumes 
for which the electron-phonon calculations were performed are summarized in Table VI. 



FIG. 11. (Left panel) Ground state energy without ZPE, including ZPE at the quasiharmonic level, and at the SSCHA level. 
The lines represent fitted curves following the recipe described in Sec. . (Right panel) The pressure derived from the fitted 
energy curves. The vertical dashed lines denote the volumes used in the SSCHA calculations. 
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FIG. 12. ZPE calculated within the quasiharmonic level for H 3 S and D 3 S. The solid lines denote the linear fit obtained with 
two data points, the volumes are those from the right panel of Fig. 11. 


TABLE V. Calculated parameters for the quadratic fit of Eq(V) and the linear fit of E v (V). 



A 0 (mev) Bo 

(meV/a.u.^) 

Co (meV/a.u. s ) 


49319 

-900.88 

3.9772 


A v (mev) B v 

(meV/a.u. 3 ) 


H 3 S quasiharmonic 

1918.5 

-12.479 


D 3 S quasiharmonic 

1404.7 

-9.0707 


H 3 S sscha 

1814.2 

-11.027 



TABLE VI. Correction to the pressure P from the vibrational energy. 


V (a.u. a ) 


P (GPa) 



no ZPE 

D 3 S quasiharmonic 

H 3 S quasiharmonic 

H 3 S SSCHA 

89.9822 

200 

206 

210 

212 

84.3035 

250 

255 

259 

261 
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